Nuclear Lattice Simulations with Chiral Effective Field Theory 

Dean Lee^, Bugra Borasoy^, and Thomas Schaefer^'^ 
^Physics Department, North Carolina State University, Raleigh, NC, 27695, USA 
^Physik Department, Technische Universitdt Miinchen D-85747, Garching, Germany 
^RIKEN-BNL Research Genter, Brookhaven National Laboratory, Upton, NY 11973, USA 

Abstract 

We study nuclear and neutron matter by combining chiral effective field theory with non- 
perturbative lattice methods. In our approach nucleons and pions are treated as point particles on 
a lattice. This allows us to probe larger volumes, lower temperatures, and greater nuclear densities 
than in lattice QCD. The low energy interactions of these particles are governed by chiral effective 
theory and operator coefficients are determined by fitting to zero temperature few-body scattering 
data. Any dependence on the lattice spacing can be understood from the renormalization group 
and absorbed by renormalizing operator coefficients. In this way we have a realistic simulation 
of many-body nuclear phenomena with no free parameters, a systematic expansion, and a clear 
theoretical connection to QCD. We present results for hot neutron matter at temperatures 20 to 
40 MeV and densities below twice nuclear matter density. 
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I. INTRODUCTION 



The nuclear many-body problem has long been recognized as one of the central questions 



,„ nuCea. ph.*. Q. TKe .adit,o„a. app.o.h to the ,na„.-bodv ptob.en, . based on the 
assumption that nucleons can be treated as non-relativistic point-particles interacting mainly 
via two-body potentials. Three-body potentials, relativistic effects, and non-nucleonic de- 
grees of freedom are assumed to give small corrections. The many-body problem is studied 
by solving the many-body Schrodinger equation. The ground-state properties of light nu- 
clei and neutron drops have been analyzed by several groups using variational methods and 
Green's function Monte Carlo QBQBHB- 

These methods have been very successful, but there are several good reasons for seeking an 
alternative approach. One reason is the desire for a theory that is more directly grounded 
in QCD. We expect this theory to explain why two-body forces are dominant, and how 
the interaction should be chosen. In addition to that, we would like to have a framework 
which allows the calculations to be systematically improved, and provides an estimate of 
the errors due to contributions that have been neglected. If we consider the interaction 
of a single nucleon with pions and external fields such a framework is provided by chiral 
perturbation theory. In a very infiuential paper Weinberg proposed to extend effective field 
theory methods to the nucleon- nucleon interaction Over the last several years effective 
ield theory methods have been applied successfully to the two and three-nucleon system 
^, Q| . Effective field theory methods have also been applied to nuclear and neutron matter, 
)ut these calculations rely on a perturbative expansion in powers of the Fermi momentum 

nlQ. 

Our aim in this work and the goal of the Nuclear Lattice Collaboration as a whole is 
to extend effective field theory methods to the nuclear many-body problem. For this purpose 
we investigate the many-body physics of low-energy nucleons and pions on the lattice. Our 
starting point is the same as that of Weinberg. We begin with the most general local La- 
grangian involving pions and low-energy nucleons consistent with translational invariance, 
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isospin symmetry, and spontaneously broken chiral symmetry. This yields an infinite set 
of possible interaction terms with increasing numbers of derivatives and/or nucleon fields. 
Degrees of freedom associated with anti-nucleons, heavier mesons such as the p, and heavier 
baryons such as the A, are integrated out. The contribution of these particles appear as 
coefficients of local terms in our pion-nucleon Lagrangian. We also integrate out nucleons 
with momenta greater than vra"^, where a is the lattice spacing. 

The operator coefficients in our effective Lagrangian are determined by fitting to 
experimentally-measured few-body nucleon scattering data at zero temperature. The de- 
pendence on the lattice spacing is described by the renormalization group and can absorbed 
by renormalizing operator coefficients. In this way we construct a realistic simulation of 
many-body nuclear phenomena with no free parameters. In our discussion we present re- 
sults for hot neutron matter at temperatures 20 to 40 MeV and densities below twice nuclear 
matter density. 

The first lattice study of nuclear matter was done by Brockmann and Frank [l^. Thev 
used a momentum lattice and analyzed the quantum hadrodynamics model of Walecka . 
Miiller, Koonin, Seki, and van Kolck lo[, were the first to look at infinite nuclear and neutron 
matter on a spatial lattice at finite density and temperature. They used an effective nucleon- 
nucleon interaction on a 4'^ lattice and found evidence for saturation in nuclear matter. The 
approach we pursue is similar in spirit to that of The main difference is the inclusion 
of pion degrees of freedom and our use of chiral effective field theory with Weinberg power 
countin g. The nuclear liquid-gas transition has also been studied using classical lattice gas 
models [iTf 



II. NOTATION 



Before describing the physics it will be helpful to first define the notation that we use 
throughout our discussion. We let n represent integer-valued lattice vectors on our 3 + 1 
dimensional space-time lattice. We use a subscripted 's' such as in to represent purely 
spatial lattice vectors. We use subscripted indices such as j for the two spin components 
of the neutron, f and [. We let be the unit lattice vector in the time direction and let 
Is = i, 2, 3 be the corresponding unit lattice vectors in the spatial directions. A summation 
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symbol such as 

E (1) 

Is 

implies a summation over values Ig — 1, 2, 3. 

We let a be the lattice spacing in the spatial direction and L be the length of the spatial 
lattice in each direction, is the lattice spacing in the temporal direction and is the 
length in the temporal direction. We let at be the ratio between lattice spacings, 

= f ■ (2) 

Throughout we use dimensionless parameters and operators, which correspond with physical 
values multiplied by the appropriate power of a. Wc use the superscript 'phys^ such as in 
m^^* to represent quantities with physical units. We use a,a^ to represent annihilation 
and creation operators for the neutron, whereas c, c* indicate the corresponding Grassmann 
variables in the path integral representation. We use the symbol : /(a^, a) : to indicate the 
normal ordering of operators in f{a^, a). We let mjv be the mass of the neutron, /i be the 
neutron chemical potential, and be the mass of the pion. 
Our conventions for Fourier transforms are 

fik)-^J2^''*-''f{n), (3) 

n 

/(^) = 7feE^"''*'"/(^)' (4) 

k 

where 

K^Cf^ko,'fh,'fk2,^ks). (5) 

We use periodic boundary conditions in the spatial directions and periodic/antiperiodic 
boundary conditions in the temporal direction for bosons/fermions. 

We let DN{k)Sij and DT^{k) be the free neutron and neutral pion propagators. For 
notational convenience the spin-conserving 5ij in the neutron propagator will from here on 

— * — * 

be imphcit. The self-energies, T,j^{k) and E^(A;), are defined by 

i^r (^) = — ^"^ffl ^ , (6) 

1 - i:^{k)D^{k) 

Di-\k) = M^''^ (7) 

1 - T.^{k)D^{k) 

where D{^^\k) and Dl^^\k) are the fully- interacting propagators. 




+ 
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FIG. 1: Chiral expansion of the two-particle irreducible kernel. 



III. NON-PERTURBATIVE EFFECTIVE FIELD THEORY 



Effective field theory provides a systematic method to compute physical observables order 
by order in the small parameter Q/A^, where is the chiral symmetry breaking scale and 
Q = {q^iTLT^, . . .). Here, g is a small external momentum and m.,^ is the mass of the pion. 
The simplest processes are those that involve only pions and external fields. In this case the 
effective field theory is perturbative. At any order in Q there are only a finite number of 
diagrams that have to be included. At lowest order these are tree diagrams with the leading 
order interaction. At higher order, diagrams with more loops or higher order terms in the 
interaction have to be taken into account. 

Weinberg showed [s] 21| that the simple diagrammatic expansion for nucleon-nucleon 



scattering is spoiled by infrared divergences. He suggested performing an expansion of 
the two-particle irreducible kernel (see Fig. Q]) and then iterating the kernel to all orders to 
produce the scattering Green's function (see Fig. Ej). It was later pointed out that a possible 
difficulty arises because at any order in Q an infinite number of diagrams is summed, and it 
is not clear that all the cutoff dependence at that order can be absorbed into counterterms 
that are present at that order 2^. This problem does indeed arise if one considers nucleon- 
nucleon scatterin g in the ^5*0 channel [2j|, but in practice the cutoff dependence appears to 
be very weak [20] j2J| . 

In this work we go one step further and consider the nuclear many-body problem. We 
expand the terms in our action order by order. 



S = S^ + Si + S2 + -- 



(8) 



At order k in the chiral expansion, we calculate observables by evaluating the functional 
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FIG. 2: The two-particle irreducible kernel is iterated to all orders to produce the nucleon-nucleon 
scattering Green's function. 

integral 

{G[N, N, vr); = - ^ ^ . (9j 

" J DNDNDn exp [—So — Si — — Sk\ 

We will refer to this approach as non-perturbative effective field theory. The interactions at 
chiral order k or less are iterated to arbitrary loop order. The functional integral is computed 
non-perturbatively by putting the pion and nucleon fields on the lattice and using Monte 
Carlo sampling. Since the number of diagrams at a given chiral order grows exponentially 
with the number of nucleons, a non-perturbative technique such as this is needed for systems 
with more than just a few nucleons. 

Computing the path integral corresponds to summing an infinite set of diagrams. As 
in the case of iterating the two-particle irreducible kernel to determine the full two-nucleon 
Green's function, it is not clear that the cutoff dependence at a given order in the low energy 
expansion can be absorbed into a finite number of coefficients in the action. In practice we 
will therefore restrict ourselves to lattice cutoffs that satisfy vra^^ < ^x- order to show 
that the effective field theory calculation is consistent we must find a window of lattice 
cutoffs such that the many-body calculation is independent of the cutoff up to terms that 
are higher order. We shall study this question numerically in our results section. 



IV. LOWEST ORDER INTERACTIONS 

Our momentum cutoff scale is 7ia~^ and we choose the lattice spacing so that 

< vra"-^ < A^. (10) 

An irreducible diagram is one that cannot be disconnected by cutting internal 
match the set of incoming or outgoing particles. In the Weinberg counting scheme 
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we estimate the chiral order of an irreducible diagram by associating one power of Q/A^ for 
each derivative interaction or exphcit factor of m,r, four powers for each loop integral, one 
inverse power for each nucleon internal line, and two inverse powers for each pion internal 
line. If V is the chiral order of an irreducible diagram, it can be shown that 

i/ = 4-f + 2/-2C+ ^ bi. (11) 

vertex i 

In (jllj) e„ is the number of external nucleons, / is the number of loops, c is the number of 
connected pieces, and bi for each vertex is 



(12) 



where j^d is the number of derivatives, j^vn.^ is the number of explicit factors of in the 
coefficient, and j^n is the number of nucleon fields. It turns out that #m,r is always an 
even number. 

We let iV represent the nucleon fields. 



iV 



proton 




T 


neutron 







(13) 



We use Ti to represent Pauli matrices acting in isospin space, and we use a to represent 
Pauli matrices acting in spin space. Pion fields are notated as tTj. We denote the pion decay 
constant as Fl^^' ^183 MeV and let 



D = 1 + T^llFl 



(14) 



The lowest order Lagrange density for low-energy pions and nucleons is given by terms with 



5^ = 



26| 



£(0) 



(Vvr, 



- D-^F-^gAN 



\2 ^2 
— TT, 



N - D-^F;^N[eijkTiTiMN 



- \Cs : NNNN ■ -\Ct : NaN ■ NaN : . 



(15) 



qa is the nucleon axial coupling, and eijk is the Levi-Civita symbol. The chemical potential 

fi controls the nucleon density and fi will be set very close to mAr. At next order we have 
terms with 5j = 1, 

^^'^ = 2^NV'N + ... (16) 
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We will include this kinetic energy term from C^^^ in our lowest-order Lagrange density so 
that we get the usual free nucleon propagator. 

In this study we limit ourselves to the interactions of neutrons and neutral pions and 
consider processes with only up to two pions. As a result we have at lowest order the terms 



-Til + (Vtto)^ + mini + a][idQ + {tun - 



(17) 



where a, a) are annihilation and creation operators for the neutron. In the Euclidean 
formalism, we have the partition function 



/ 



DnDNDNexp 



DnDNDNexp / d'^x Ce 



(18) 



where 



^0 + (Vtto)^ + mini - a][do + {m^ - ^)]aj 



+ ^a\aijaj ■ Vtto — Ca|a-|-ajaj^. 



(19) 



We will use Xq to represent the Euclidean temporal coordinate, rather than switching from 

Xq to Xi. 

V. FREE NEUTRON 



In the simplest discretization, the Euclidean lattice action for free neutrons has the form 



-^^T'^ = XI + 6) + (-1 + (mjv - ij)at + 6/i) c*{n)ci{n)\ 

n,i 

-hJ2 [cmciin + is) + c*{n)ci{n - I)] (20) 

where 

h = ^- (21) 

However we want a discretization that will minimize the dependence on at so that fewer 
lattice steps in the temporal direction can be used and results for different at can be directly 
compared. 



Let us review the conversion from the operator formahsm to path integrals. The free 
neutron lattice Hamiltonian is 



Hnn = Yl [("^^ ~ ^ + 7^)4{ns)ai{ns) 

ns,i 



ns,ls,i 

We want to convert the partition function for free neutrons, 

Z = Tt [exp(-/3if^^)] = Tr [exp^-atHpf^) exp(-atifjviv)--- exp(-atif^Ar)] 
into the form 



(22) 



(23) 



Using the identity 



we can write 



Z = J DcDc* exp [-5/ 
exjp[alXijaj] =: exp [a J (e^ — l)ijaj] :, 
exp(-Q;ji/7VAf) =: exp{-hf^^{a' , a)) : +0{h), 



(24) 
(25) 
(26) 



where 



_ /,e-{(-iv-M)"*) J2 \al{n,)a,in, + l) + alins)a^in,-l) 



(27) 



Introducing the extra e ^^"^^ factor multiplying h is not well motivated at this stage, 
but it insures that the neutron chemical potential is coupled to an exactly conserved neutron 



number operator 



28| 



. We now use the correspondence 



2210 



Tr [: /„_i(a^a) : ... : fi{a\a) :: /o(a^a) :] 



= dCn-idcl_^...dcodc*Qexp[ c*{cj - Cj+i)] JJ fjic*,Cj) 

i=0,...,n-l i=0,...,n-l 

with c„ = —Co- We can now convert the partition function to the path integral form in 
with 



(28) 



Snn = Yl + 0) - e 

ns,i 



l^c* (n)Q(n + /,,)+ c*{n)ci{n - Is 



(29) 



nisi 
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This lattice action has temporal discretization errors of 0{h), whereas the action in ()2Up has 
errors of 0{at). Since /i is a small parameter this is an improvement and so the dependence 
on at has been significantly reduced. 

It is conventional to define a new normalization for q, 

-(m]v-/i)at 



Then 



where 



= c^e-^"^'^-^^"*. (30) 
Z = e-2{™iv-M)/3i^ J i)c'Dc* exp [Sm] (31) 

Snn = E [e(™--^)"*c*(n)c:(n + 0) - e-''c*{ny,{n)] 

n.i 

-hJ2 [cKn)c^{n + Is) + c*{nyi{n - Is)] • (32) 

n,la,i 

We observe that the neutron chemical potential is in fact coupled to an exactly conserved 
neutron number operator since it appears in the same manner as a lattice gauge connection 
in the temporal direction. Comparing the two actions ()20|) and ()32|) we can summarize the 
difference as follows. If we write 



where 



h ^ l^c* (n)ci(n + Is) + c*{n)ci{n - 4) 

n'ls.i 



(33) 



then 



X = 6h, 

Sm = E [e('"--^)-c*(n)c:(n + 0)] - J] [c*(n)e-^c:(rl)] 

n,i n,i 

- h [cm^in + Q + c:{n)4{n - I)] . 



(34) 



(35) 
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In momentum space we have 



^NN — 5^ ^)c^(^) 
k,i 



We now have the free neutron correlation function on the lattice, 



jDc'DcXin)cm^M-SNN] _ 1 



/ Dc'Dc* exp [Snj^] LtL^ 
(no sum over i) where the free neutron propagator is 
DN{k) ^ 



VI. NEUTRON CORRELATION FUNCTIONS 



(36) 



(37) 



(38) 



At nonzero time step there are some subtleties going from the correlation functions in 
the operator formalism to correlation functions in the path integral formalism. We have 



Z = Tr [exp(— a^i?)... exp{—atH) exp{—atH)] 
= Tr [: /L(-l(a^ a) : • • • : fi{a\ a) :: /o(a^ a) :] 

where the total number of time steps is Lt and each of the /^ 's are the same, 

: /j(a^ a) :— exp(—atH) + 0{h). 

Suppose we wish to calculate 

Tr [uLt-i{a)vLt-i{a^) : ht-i ■ ■ • • ui{a)vi{a') : fi : uo{a)vo{a') : fo :] 
This can be rewritten as 



(39) 



(40) 



(41) 



Tr [: VLt-i{a^)fLt-iULt-2{a) vi{a^)fiUo{a) :: vo{a^)foULt-i{a) :] 

In the path integral formalism this is equivalent to 

J DcDc*F{c, c*)exp[-,S] 

where 

F{C,C*) = VLt-lic%_i)uLt-2{cLt-l) ' ' ' Vo{cl)uLt-l{cQ) ■ 

We note that Vj is a function of c* whereas Uj is a function of c*,^^ 
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(42) 



(43) 



(44) 



VII. FREE NEUTRAL PION 



The lattice action for a free neutral pion is 

Snn = ((^ + 3)a;t + at^) ^ 7r(n)7r(n) - ^ [ei7r(n)7r(n + /) 

n 

where 



n,l 



(eo, ei, 62, 63) = (q;^ \ a^, a^). 
In momentum space the action is 



and so 



J Dmi{n)ii{0) exp [—8^, 



J Dfrexp [—Sti 
where free neutral pion propagator is 

dM = 



-y 

.1.3 



(45) 



(46) 



(47) 



(48) 



(49) 



In this first exploratory study we are not concerned with the issue of exact chiral symmetry 
on the lattice and therefore will neglect the Haar measure. This aspect of exact chiral 
symmetry will be investigated in a future study along with the inclusion of charged nucleons 
and pions. 



VIII. NEUTRAL PION-NEUTRON COUPLING 

In the continuum the pion-nucleon coupling makes a contribution to the integrand of the 
partition function of the form 



exp 



(50) 



where a — ((71,(72,(73) are Pauli matrices for spin, = (ti,T2,T3) are Pauli matrices for 
isospin. 



(51) 
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and FP'^y^ ^183 MeV is the pion decay constant, 



We keep only the term involving the neutral pion and neutron, 



exp 



^ I d'^x ( c*aijCj ■ Vtto 



The simplest lattice discretization of this interaction term is 



where 



exp [S^nn] 



(52) 



(53) 



(54) 



and 



n 

n 
n 

A^='=7ro(ri) = 7ro(n + Z) - 7ro(ri - I). 



(55) 



(56) 



Wc now use a temporally-improved discretization. We can write the simple lattice action 
for the free neutron with pion-neutron coupling as 

n,i 

+ E + - l^)oit]5ij + Xij{n)) Cj{n) 

- ^ E (^)ci(^ + L) + c*{n)ci{n - Is) (57) 



n,ls,i 



with Xij{n) given by the matrix 



■|^(Af + zA±)7ro(n) 



9Am 



(Af -^A±)7ro(n) 



(58) 
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Then our temporally-improved action is 



n] 



- \c*in)Ciin + ls) + c*in)Ciin - Is 

nis.i 



(59) 



IX. NEUTRON CONTACT TERM 



The neutron contact term has the form 



Hnnnn = C^a\{n)a^{n)a\{n)ai{n). 



(60) 



We can rewrite the contribution at lattice site n to the partition function using a discrete 
Hubbard-Stratonovich transformation Isill. For C < 0, 



exp(— Ca^aj {n)a^ (ri)aj {n)ai (n)) 



i J2 exp 

s{n)=±l 



■ Cat 
. 2 



+ As(n))(a|(n)a|(n) + a|(n)aj(r?) — 1] 



where 



Since 



we can write 



cosh A = exp(— ^Y^-). 



+ e"^ = 2exp(-- 



2 ^' 



A = ln 



exp 



+ v/exp(-Cat) - 1 



The simplest lattice discretization gives a contribution to action, 



(61) 

(62) 
(63) 
(64) 



where 



and 



C/3 r 3 I c _|_ Qsimple 



Sss = As(n) 



(65) 
(66) 

(67) 



15 



However this actually gives a result that is inconsistent with the Hamiltonian operator 
form in ()6ip in limit at 0. The problem is somewhat subtle and has not been given 
much attention in the literature. The point is that operator ordering at O(A^) cannot be 
ignored since ~ 0{at). We will deal with this in the same way that we constructed the 
temporally-improved action for the free neutron and the pion-neutron coupling. We write 

rfsimple , Qsimple , Qsimple 
^NN '-'ttNN ^sNN 

= ^c*{n)ci{n + 6) 

+ ^c*{n) ([-1 + {mN - lJ^)at]Sij + Xij{n)) Cj{n) 



- \c*{n)ci{n + Is) + c*{n)ci{n - Ig) 



(68) 



with Xij{n) equal to 



6h - 2^Atvro(n) + ^ + A 



(Af + zAf ) 7ro( 



s n 



9Aat 
2F^ 



{Af - lAf) vTol 



n) 



n] 



(69) 



then 



S 



NN+ttNN+sNN 



E 



c*(n)c:(n + 0) 
J2cm (e^^("))^^.c;.(n) 



h ^ ^c*{n)ci{n + 4) + c*{n)ci{n - 4) 



(70) 



X. ONE-LOOP NEUTRON SELF-ENERGY 



In the next few sections we calculate several lattice Feynman diagrams. These calcula- 
tions will serve as a check that our non-perturbative simulation is functioning properly in 
the small coupling limit ^ 0, C — > 0. It will also give us a reference point to measure 
how non-perturbative the interactions are at physical values for qa and C and for various 
densities. 

At 0{g\) we have a contribution to the neutron self-energy due to a neutral pion and 
neutron intermediate state as shown in Fig. |21 If we write the pion-interaction interaction 
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k+p, j 



FIG. 3: Neutron self-energy due to a neutral pion and neutron intermediate states. 



k, i 




k, i 



FIG. 4: Neutron self-energy due to the ttttNN interaction in the temporally-improved action. 



term to order 0{gA) in momentum space, we have a contribution to the action of the form 



r5Z - p)d^ik) - c\{-k - p)c!fi) 7ro(p) sin(p^3) 



k,p 



+ W^Y. m-k-p)~dfi) 7ro(p)[sin(p.i)-zsin(p.2 



+ 



k,p 

fwl-^lI Y1 [^U~^ ~ P)c\{k) TToip) [sin(p,i) + i sin(p,2)] • 

k,p 



(71) 



Then the diagram in Fig. El leads to a contribution to the self-energy that goes as 

Ef\k) = ^^f^Y.DMDN$ + p) [sin2(p,0 +sin2(p,2) + sin^(p*3)] • (72) 

p 

Our temporally-improved action has a ttttNN interaction of the form 



S = ...- ^e"^^ cm4in)iAtMn)f 



ftls.i 



(73) 



which gives rise to the diagram in Fig. 0] We get an additional contribution to the self- 
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k, i 




k, i 



energy, 



where 



FIG. 5: Neutron self-energy due to the contact interaction. 



(74) 



(75) 



At 0{C) we have the one- loop diagram shown in Fig. El If the vertex is located at lattice 
site n we can isolate the relevant lowest-order interaction in the path integral starting from 



^{Cat/2+Xs{n)) 



s{n)=±l 

Expanding the exponential we get 



+cj(n)c'^(n) 



5 E ^ 

s{n)=±l 



^g-i2h (g-(Cat/2+A.(n)) _ l) " c| (n)^ (n)cj (n)c'^ (n) 



We find that 



(76) 



(77) 



s{n)=±l 



{Cat/2+\s{n)) ^^-{Cat/2+Xs{n)) _ ^ g 



and 



s(n)=±l 

So to lowest order we can write the interaction as 



(78) 



(79) 



-12h (^^-Ca, _ c*(n)c;(n)cj(n)c'^(n). 



(80) 
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k+p, j 

FIG. 6: Pion self-energy due to neutron-neutron hole intermediate states. 

Therefore the contribution to the self-energy is 

e'^P (k) = e"i2'^(e"^"* - (81) 

where 

/Dc-Dc*cKO)c^(0)exp [-%^] 

Y = — (no sum over t) 

J Dc'Dc* exp [-D^TvJ 

^E^^(^)- (82) 



LtL^ . 

k 



This self-energy term is proportional to Y ^ which in turn is proportional to the neutron 
density with an 0{at) time discretization correction. 

XI. ONE-LOOP PION SELF-ENERGY 

At 0((7^) we have a contribution to the pion self-energy due to a neutron and neutron- hole 
intermediate state as shown in Fig. IHl The contribution to the self-energy is 

T.^"^\k) = - '"g:;"' DN{p)DN{k + p) [sin2(A;,i) + sin2(A;,2) + sin2(A;,3)] , (83) 

p 

Our temporally-improved action has a term of the form 

5 = ■ ■ ■ - '^e-'^ ^ c*(n)c:(n)(Ajvro(n))^ (84) 

This leads to the diagram in Fig. [7| and gives an additional contribution 

T.f^\k) = ... + ^e-6^y [sin2(fc,i) + sin2(fc,2) + sin2(A;,3)] • (85) 
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p, i 

C) 



FIG. 7: Pion self-energy due to the ttttNN interaction in the temporany-improved action. 




k+p, j 

FIG. 8: Two loop connected bubble diagram at 0{g\). 



XII. TWO-LOOP AVERAGE ENERGY 

We will calculate the shift in the average energy by computing 

The logarithm of the full partition function is the sum of the connected diagrams. At 0{g\) 
we get a contribution from the connected bubble diagram shown in Fig. |H1 

The amplitude for this bubble can be obtained in a straightforward manner from either 
(|72|) or (jHl: 



F^L'iLt 5Z DNip)DNik + p)D^{k) [sin^(A;*i) + sin^iKi) + sin^(/c*; 

p,k 

In our temporally-improved action the term 

5 = ... - ^e-^'^ Yl cmc'MiAtMn))'. 



(87) 



(88) 



nls.i 



produces the diagram shown in Fig. 121 The amplitude for this process can be computed 
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FIG. 9: Two loop connected bubble due to a mrNN interaction in the temporally-improved action 
at 0{g\). 




FIG. 10: Two loop connected bubble at 0{C). 



from either (j74j) or (j85p and is 

4#e-^'^F DAk) [sm\K,) + sin2(A;,2) + sin2(A;*3)] = ^^^^e-^^XF. (89) 



9a^ 



At 0{C) we have the connected bubble diagram shown in Fig. ^| If the vertex is located 
at lattice site n the lowest order interaction is 



-i2h (^^-ca, _ c^nmn)clin)c\{n). 



(90) 



Summing over all sites, we find that the connected bubble in Fig. ^Ihas the amplitude 



(91) 
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XIII. AVERAGE ENERGY FROM SIMULATIONS 



The average energy can be computed by taking —-^InZ and then adding /lA to the 
result, where A is the average number of nucleons. The partition function is given by 

Z{P) oc e-2(---'')/^^%-^' / DttDsDc'Dc* exp [-%iv+.iviv+.iviv - - Sss] ■ (92) 

It is convenient to define a new partition function Z' with normahzation, 

Z'iP) = J DnDsDc'Dc* exp [--SjvAr+.iviv+.iViv - - Sss] ■ (93) 

Z' is what we actually compute in the simulation. Then 



d_ 
"dp 



\nZ 



d_ 
"dp 



In Z' + 2{mN - l^)L' 



(94) 



After computing — ^ In Z we subtract out the value for —^\nZ at the same P but zero 
neutron density. In our theory we have not included pion self-interactions. Therefore in 
the absence of neutrons we can calculate —-^InZ for free pions exactly. If we later decide 
to include pion self-interactions then a separate simulation with only pions will be needed 
for this calculation. 

At zero neutron density and free neutral pions, the lattice path integral gives 

-1/2 



Z' oc 



det 



(95) 



where (>S'^7r)jj ai'e the coefficients of the quadratic form in the pion action S^tt- We find 



^ 1 



[det ((5, 



-1/2 



^dp 



+ 3)at + a^^ - ^ei cos{Ki) 



(96) 



This quantity is the pion contribution to —-^InZ'. This is not exactly what one might 
define as the pion contribution to the energy, though it is closely related. To calculate the 
pion contribution to the energy, one also needs to add 

d_ 
dp 

which arises from an additional factor of 



WLtP-\ 



(97) 



(98) 
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in the partition function. This factor is due to the conjugate momenta integrations going 
from the Hamiltonian formahsm to the path integral formahsm, 



dp (g„+i| exp 



AT 

2 d^q 



\p) {p\ Qn) oc J dpexp [-^p^ + ip{qn+i - Qn 



After we calculate —-^InZ, we can compute the average energy per nucleon using 

A A dp ^ 
In the simulations Lt is kept fixed, and at is varied in order to compute 



(99) 



(100) 



XIV. WEAK COUPLING RESULTS 

In this section we check that the results of our numerical simulations at weak coupling 
agree with our results from perturbation theory. For the neutrons we define the temporal 
two-point correlation functions, 

(a^a^) = Z-^Tr [exp(-(/5 - ntat)H)a^{nt, 0, 0, 0) exp{-ntatH)a\{0, 0, 0, 0)1 (101) 

and spatially separated correlation functions in the x direction 



fa a^' 



) = Z-^Tr exp{-(3H)a^{0, n„ 0, 0)a{(0, 0, 0, 0) 



(102) 



The results are exactly the same in the y and z directions. Similarly for the neutral pion, 
we define the temporal and spatial correlation functions 



{n^n) = Z-^Tr [exp(-(/5 - ntat)H)iT{nt, 0, 0, 0) exp{-ntatH)7r{0, 0, 0, 0)] , (103) 



(ttJ^tt) = Z-^Tr[exp(-^i7)7r(0,n„0,0)7r(0, 0,0,0)]. (104) 

At weak coupling we compare the temporal and spatial correlation functions for the 
neutron and pion as well as the energy per neutron, E/A. We use the parameters = 150 
MeV, p = 2.0, L = 3, Lt = 3, mf^' = 939 MeV, and = 135 MeV, /i = m^v - 0.1. We 

first take gA = ^ and C = —0.135. In Tables la and lb we show the results for the free 
neutron correlation functions, the one-loop results using the self-energy correction given in 
(|HT|). and the results of our Monte Carlo simulation. 
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Table la: 



;a at) for c/a = 0, C 



-0.135 



nt 





1 


2 


Free 


0.7568 


0.5027 


0.3444 


One-loop 


0.7453 


0.5059 


0.3537 


Simulation 


0.7447(2) 


0.5057(3) 


0.3537(3) 



Table lb: (a^a^) for c/^ = 0, C = -0.135 



Us 


1 


Free 


-0.03903 


One-loop 


-0.03940 


Simulation 


-0.03936(2) 



In Table Ic we show the free neutron value for the energy per neutron, the two-loop corrected 
value using (jHH), and the result of the simulation. 

Table Ic: E/A for = 0, C = -0.135 



Free 


6.665 


One-loop 


6.653 


Simulation 


6.652(1) 



There is no correction to the free pion correlation function when qa = 0. We see that all 
the simulation results match the loop calculations for Qa = and small C. 

For C = and small qa there are two sets of diagrams which we would like to separately 
compare with simulation results. The temporally-improved neutron action has the form 



S 



NN+ttNN 



^e(™^-^)"*c*(n)c^(n + 0) 



X{n) 



^1,1,0 



h ^ l^c* (n)Q(n + is) + c*{n)ci{n - is) 



(105) 



where Xij{n) 



6h 



2^Af7ro(n) 



9-4 



2F^ 

^(A± + zAf)7rol 



n] 



(Af-.Af)vro( 
6/i + %^Af7ro(n) 



n] 



(106) 



The one- loop corrected neutron correlator gets a contribution from both ()72|) and ()74|) : the 
one-loop corrected pion correlator uses (|HHj) and (|H^: and the one-loop corrected energy per 
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neutron has terms ()87p and ()89p . The comparisons with simulation results for qa = 0.750, 
C = 0, are shown in Tables 2a-e and are labelled by 'exp', which stands for the exponential 
form used in the temporally-improved action. 

We will also remove the temporally-improved diagrams which gave us the contributions 
flTij) . (|H^. and (jSni)- We do this by replacing the term in the action 



by 



where Mjj(n) is 



exp(— 6/i) 



n,i,3 



X{n) 



-1 - S^A^TToin) 



9Aat 



f (A± - ZA±) TTol 



(Af + ^Af ) 7ro(n) -1 + ^A±7ro(rr) 



n] 



(107) 



(108) 



(109) 



The comparisons with simulation results for this linearized action are also shown in Tables 
2a-e and are labelled by 'lin'. 



Table 2a: {a^a)) for qa = 0.750, C = 



nt 





1 


2 


Free 


0.7568 


0.5027 


0.3444 


One-loop (lin) 


0.7586 


0.4978 


0.3400 


Simulation (lin) 


0.7585(1) 


0.4974(1) 


0.3399(1) 


One-loop (exp) 


0.7496 


0.5005 


0.3475 


Simulation (exp) 


0.7494(1) 


0.5000(1) 


0.3472(1) 



Table 2b: (a^a^) for qa = 0.750, C = 



Us 


1 


Free 


-0.03903 


One-loop (lin) 


-0.03859 


Simulation (lin) 


-0.03856(1) 


One-loop (exp) 


-0.03890 


Simulation (exp) 


-0.03886(1) 
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Table 2c: 



vr vr 



) for QA = 0.750, C = 



nt 





1 


Free 


0.1764 


0.0615 


One-loop (lin) 


0.1780 


0.0644 


Simulation (lin) 


0.1780(3) 


0.0643(3) 


One-loop (exp) 


0.1810 


0.0660 


Simulation (exp) 


0.1809(2) 


0.0659(2) 



Table 2d: 



) for gA = 0.750, C = 



ris 


1 


Free 


0.0364 


One-loop (lin) 


0.0366 


Simulation (lin) 


0.0364(2) 


One-loop (exp) 


0.0367 


Simulation (exp) 


0.0365(2) 



Table 2e: E/A for qa = 0.750, C = 



Free 


6.665 


One-loop (lin) 


6.673 


Simulation (lin) 


6.672(1) 


One-loop (exp) 


6.640 


Simulation (exp) 


6.638(1) 



We see that all the simulation results match the loop calculations for C = and small qa- 



XV. RENORMALIZATION OF COEFFICIENTS 



We now discuss the renormalization of operator coefficients in our lowest order effective 
Lagrangian. At zero temperature and yU < m-Ar, the pion self-energy vanishes since there are 
no neutron holes. Thus there is no renormalization for the pion wavefunction and mass. 

In the Weinberg counting scheme, the neutron self-energy at zero temperature and /i < 
m^v gets a contribution from diagrams such as the one shown in Fig. |21 This is the lowest 
order diagram and comes at chiral order v = Since we require cutoff independence, the 
counterterm diagrams must also be at order z/ = 3. Since this is a small correction, we 
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FIG. 11: One- loop correction to the one-particle irreducible ttNN vertex. 



will ignore wavefunction and kinetic energy renormalization for the neutron in the present 
study. Although the mass counterterm is also small, we will take some extra care with 
this one since we are interested in precise measurements of the energy per neutron. In the 
non-relativistic formalism the mass counterterm can be regarded as a shift in the definition 
of the chemical potential. Its purpose is to eliminate cutoff dependence in loop diagrams, 
but we will also use it to absorb residual effects due to the finite temporal spacing. 



We will refer to the mass counterterm as Am at. In the limit of zero neutron density, the 
neutrons behave as free particles. We can therefore calculate Am^r in that limit. From 
a theoretical point of view it would be nice to measure the average energy at both zero 
density and zero temperature. Computationally, however, it is more practical to make the 
measurement at non-zero temperature. 

At zero temperature and yU < mjv, the one- loop contribution to the one-particle irre- 
ducible 71 NN vertex is shown in Fig. ^2 This process is also at chiral order z/ = 3. Since 
this is a small correction we will also ignore renormalization of the ttNN coefficient and set 
its value to equal the physically measured value. 



At zero temperature and fi < rriN, the lowest order contribution to the A^A^ scattering 
Green's function is due to iterating the lowest order two-particle irreducible diagrams. The 
lowest order two-particle irreducible diagrams are shown in Figs. El and El When there is 
a bound state or the scattering length is very large compared to other relevant scales, the 
two-particle irreducible kernel must be iterated and summed to all orders. The result is 
that the A^A^ scattering Green's function has a cutoff dependence that cannot be regarded 



a* = ^ > 0. 



(110) 



9A 



1.25. 
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FIG. 12: Two-particle irreducible one pion exchange diagram. 



FIG. 13: Two-particle irreducible contact diagram. 

as a small correction. We will therefore need a non-perturbative calculation of the NNNN 
contact interaction counterterm. We will use the Schrodinger equation on the lattice to 
deal with this problem, and we describe the procedure in the next section. 

XVI. LATTICE SCHRODINGER EQUATION AND PHASE SHIFTS 

We adjust C, the coefficient of the NNNN contact interaction, so that the A^A^ s-wave 
scattering length matches the experimental value (see for example 132] ). In order to calculate 
the phase shifts, we will solve the lattice Schrodinger equation for the two-neutron system 
and observe the asymptotic form of the scattering wavefunctions. The ffist step will be to 
construct the potential between two neutrons. 

Let |0) be the free non-interacting vacuum. The two- neutron state with zero total spatial 
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momentum, zero total intrinsic spin, and spatial separation can be constructed as 

1^^) = TT^ff + Ths)a\{ms) - a|(ng + rhs)a\{ms) |0) . (112) 

m 

We let V2N be the lowest-order potential in the Weinberg counting scheme between two- 
neutron states with zero total spatial momentum, zero total intrinsic spin, and spatial sep- 
aration fis. Using the fact that 

^^^1 • 4''i^s)^f'j''(^j'i^l)4'KWij'(^j'i^'s) ■ \^s) = {ns \ns) , (113) 

we have 

T/ (-^ 3a ^ e-^"--^^-E^.sin^(fcs.)^. , 
^^"(«.) - ^ E 4,3-E.cos(...), + 

After obtaining l^2iv we can construct a matrix representation for the Hamiltonian in the 
two-neutron sector and solve the time-independent lattice Schrodinger equation. At this 
stage one could implement Liischer's formula for the measuring phase shifts in a cubical 
periodic bo. BH. Howeve. i„ ou. case we can cons.uet the ei,envec.o.s explicitly 
using Lanczos iteration, and so we find it more straightforward and accurate to read the 
phaseshifts directly from the asymptotic forms of the s-wave scattering states. Since we 
are working in a periodic box it is important to measure the phase shifts far away from the 
center of the potential and all its translations in the periodic box, as shown in Fig. ITU 

Before using this technique for our actual neutron system, we first test our technique for 
hard sphere scattering where the exact result for the s-wave phase shifts is well known. If 
the spheres have radii r/2 (therefore the centers of the spheres are separated by r) then the 
s-wave phase shift has the form 

So = -kr, (115) 

where k is the momentum. The momentum k for a given scattering state can be determined 
from the energy and the free particle non-relativistic dispersion relation. In Table 3 we show 
results for the s-wave phase shift at inverse lattice spacing = 150 MeV, with particle 
masses set at mjy and a lattice volume of 20^. We use several radii r and show comparisons 
with the exact continuum result —kr. The results suggest that our lattice Schrodinger 
technique seems to be functioning properly. 
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/ r' 


- 


. 1 


~ sin(foc+5jj(i))/x ~ sin(fcc'+8j|(^))/x' |^ 









FIG. 14: Measuring s-wave phase shifts from the asymptotic form of scattering states in a periodic 
box of length L. The center of the potential is at x = and x = L 



Table 3: Hard sphere scattering phase shifts 



r 


k 


So 


—kr 


4.5 


0.117 


-0.56 


-0.53 


5.5 


0.147 


-0.85 


-0.81 


6.5 


0.179 


-1.20 


-1.16 



We now use the lattice Schrodinger technique to tune the couphng C to reproduce the 
large scattering length that is observed in nature. In Table 4 we show the best fit values 
for C'^^^'^ for several different lattice spacings. 

Table 4: C^^^"^ for different lattice spacings 





^phys 


150 MeV 


-4.0 ■ 10-5 MeV-2 


200 MeV 


-3.4 ■ 10-5 MeV-2 


250 MeV 


-3.1 ■ 10-5 MeV-2 


300 MeV 


-2.9 • 10-5 MeV-2 



We can compare this with the pionless case where the only interaction is the contact interac- 
tion. In this case we need only sum bubble diagrams and that gives the relation C^^^^^ oc a. 
The pionless calculation on the lattice has been discussed in j35| . 
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XVII. ZONE DETERMINANT METHOD 



It can be shown that fermions at inverse temperature (5 with spatial hopping parameter 
h' have a locahzation length Q of 

/ ~ vW- (116) 

This idea was used in js^] to generate an algorithm called the zone determinant method to 
speed up the calculation of determinants using LU decomposition in nuclear lattice simula- 
tions. 

The technique is relatively simple to describe. Let M be the neutron matrix, in general 
an n X n complex matrix. We partition the lattice spatially into separate zones such that the 
length of each zone is larger than the localization length I. Since most neutron worldlines 
do not cross the zone boundaries, they would not be affected if we set the zone boundary 
hopping terms to zero. Hence we anticipate that the determinant of M can be approximated 
by the product of the submatrix determinants for each spatial zone. 

Let us partition the lattice into spatial zones labelled by index j. Let {Pj} be a complete 
set of matrix projection operators that project onto the lattice sites within spatial zone j. 
We can write 

M = Y^ PjMPi = Mo + Me, (117) 

id 

where 

M^ = Y,P^MPi, (118) 

i 

If the zones can be sorted into even and odd sets so that 

PjMPi = (120) 

whenever i is even and j is odd or vice- versa, then we say that the zone partitioning is 
bipartite. We now have 

det(M) = det(Mo) det(l + Mq^Me) 

= det(Mo) exp(Tr(log(l + Mq^Me))). (121) 
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Using an expansion for the logarithm, we have 



det(M) = det(Mo) exp [ ^ ^ ^ TT{{MQ^MEy) ) . (122) 



Let us define 



Am = det(Mo) exp ^-^—Tr{{M^^MEr)j ■ (123) 

Let \k{MQ^ Me) be the eigenvalues of Mq^Me and R be the spectral radius, 

R= max (|Afc(Mo"^M£;)|). (124) 

k=l,...,n 



(125) 



It has been shown M that for < 1 



|det(M)-A^| ^ 



where 

c= -nlog(l-i?). (126) 

The spectral radius R determines the convergence of our expansion. R can be reduced by 
increasing the size of the spatial zone relative to the localization length /. In the special 
case where the zone partitioning is bipartite, we note that for any odd p, 

1:i{{M^^MeY) = 0. (127) 

In that case 

^2m+l = ^2m- (128) 

In the simulations presented in this article we use the zone determinant method to calcu- 
late neutron matrix determinants. We use the second order approximation A2 with zones of 
the smallest possible size, a single spatial point ([1,1,1] in the notation of j..37]). An estimate 
of the approximation error is discussed along with each measurement in the results section. 



XVIII. RESULTS 

We have generated simulation results for a^^ = 150 MeV; at = ^ = 1.0; temperatures 
j-phys = 25.0 MeV and 37.5 MeV; and lattice sizes 3^ 4^, and 5^ Half-filling at this lattice 
spacing occurs at 

— = 2.64, (129) 
Pn 
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FIG. 15: Energy per nucleon in MeV for temperatures 25 MeV and 37.5 MeV and different lattice 
volumes. 

where p^^* is the normal nuclear density of about 0.17 nucleons per fm^. The calculations 
were performed using the zone determinant method using the second order approximation 
A2 with zones consisting of a single spatial point. By calculating the exact determinants of 
some generated matrix configurations we estimate the systemic error for the zone expansion 
to be about < 0.1% for Tp^^' = 37.5 MeV and < 0.5% for TP'^y = 25.0 MeV. 

We have dealt with the complex action by computing the phase as an observable, 

For the various simulations presented here we found an average phase of about 

g- Re5[n]g— ilm5[n] 

0.95 - 1.00, (131) 

and so this did not present a significant computational problem. 

In Fig. El we show the energy per neutron as a function of neutron density. Our results 
indicate a rather fiat function for the energy per neutron as a function of density. For 
comparison we show in Fig. El the energy per neutron for the free neutron for temperatures 
rpphys ^ 7 ^5 25.0, and 37.5 MeV. 

In Figs. El and El we compare with results for free neutrons on the lattice and loop 
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Energy per neutron for a' ^ 1 50 MeV, L ^ 4 
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FIG. 16: Energy per nucleon in MeV for temperatures 7.5, 15, 25 and 37.5 MeV and lattice volume 
43. 



calculations for the physical values of qa and C. As is easily seen, the loop calculations are 
not very close to the non-perturbative simulation results for the densities shown. 

The flatness of our energy per neutron curves at these temperatures are intriguing and 
hopefully will be checked by others in the near future. Our results appear to be consistent 
with the neutron matter results of at the same temperatures. Variational calculations 
js^ and recent quantum Monte Carlo results from 2] observe a gradually flattening of the 
energy per neutron curve with increasing density, though not as flat as the results we see. 
The calculations in 1^ however were performed at zero temperature and at lower densities. 

We implemented a temporally-improved action in order to remove as much as possible 
the dependence on = ^, the ratio of the temporal lattice spacing to the spatial lattice 
spacing. In Fig. UHlwe show the dependence on at for at = 1.0, 0.667, 0.500. We see that 
the dependence on at is minimal. 

We now look at how the energy per neutron changes as the interaction strength is varied. 
According to Table 4, the physical value for CP'^y at a"^ = 150 MeV is -4.0 ■ 10"^ MeV~^ 
If we however take the coupling to be 50%, 100%, and 150% of the physical value we find 
the results shown in Fig. [201 

In Figs. 1^ and 1221 we show density versus chemical potential and comparisons with the 
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Energy per neutron for = 150 MeV, L = 3 at 7*"^'^ = 37.5 MeV 
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FIG. 17: Energy per neutron for temperature 37.5 MeV and comparisons with the free neutron 
result on the lattice and loop calculations. 



free neutron and loop calculations. We can see again that loop calculations are not close to 
the simulation results. For fixed chemical potential, the simulation density is higher than 
the loop-calculated density, which is in turn higher than the free neutron density. 

If our effective field theory formalism is valid, we should be able to reproduce the same 
results for different lattice spacings. There are however practical computational constraints 
on a. The determinant zone expansion for fixed zone size and expansion order will break 
down if the lattice spacing is too small. For a'^ = 200 MeV we estimate that the second 
order approximation A2 with zones consisting of a single spatial point produces errors of size 
roughly 3%. In Fig. 1221 we compare the energy per neutron as measured for = 150 MeV 
and a^^ = 200 MeV at temperature 37.5 MeV. The results are in rather good agreement. 
In the future we hope to use larger zone sizes and do simulations at lattice spacings up to 
a-^ = 300 MeV. 

XIX. SUMMARY AND FUTURE DIRECTIONS 

We have introduced a new approach to the study of nuclear and neutron matter which 
combines chiral effective field theory and lattice methods. Nucleons and pions are treated on 
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Energy per neutron for a"^ = 1 50 MeV, L = 3 at T^'^''^ = 25.0 MeV 
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FIG. 18: Energy per neutron for temperature 25.0 MeV and comparisons with the free neutron 
result on the lattice and loop calculations. 
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Energy per neutron for a'^ - 1 50 MeV, L = 3 
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FIG. 19: Dependence on at for temperature 37.5 MeV. 



the lattice as point particles and wc arc able to probe larger volumes, lower temperatures, and 
greater nuclear densities than in lattice QCD. The low energy interactions of these particles 
are governed by chiral effective field theory and operator coefficients are determined by 
fitting to nucleon scattering data. Any dependence on the lattice spacing can be absorbed 
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Energy per neutron for a"^ = 1 50 MeV, L = 3 at T^'^''^ = 37.5 MeV 
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FIG. 20: Dependence on CP^^y^ for temperature 37.5 MeV. 



Density vs. chemical potential for a'^ = 1 50 MeV, L = 3 at 7*''^'^ = 37.5 MeV 
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FIG. 21: Density versus chemical potential for temperature 37.5 MeV. 



by the rcnormalization of operator coefficients. In this way we have a realistic simulation 
of many-body nuclear phenomena with no free parameters, a systematic expansion, and a 
clear theoretical connection to QCD. We have presented results for the energy per neutron 
for hot neutron matter at temperatures 20 to 40 MeV and densities below twice nuclear 
matter density. 
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Density vs. chemical potentiai for a"^ = 1 50 N/leV, L = 3 at 7*''^^^ = 25.0 H/leV 
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FIG. 22: Density versus chemical potential for temperature 25.0 MeV. 



Energy per neutron for 7*''''^ = 37.5 MeV 
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FIG. 23: Comparison of the energy per neutron results for a ^ = 150 MeV and a ^ = 200 MeV at 
temperature 37.5 MeV. 



In conjunction with other members of the Nuclear Lattice Collaboration, we plan several 
extensions, generalizations, and improvements upon this work. In the course of producing 
data for this article, we have also generated a large amount of data for the neutral pion, 
neutron, and pair density correlation functions. This data will be analyzed and presented in 
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a forthcoming article. In the near future we also plan to study the pionless version of the same 
neutron system. There has been a recent mean- field discussion of this model . Without 
pions, the phase will be completely eliminated from the matrix determinant. This is due to 
the fact that the matrix is purely real, and in the Hubbard- Stratonovich formalism the up 
and down spins appear in a way that the matrix determinant is the square of a real number. 
We also plan to extend our studies to include neutrons, protons, neutral and charged pions, 
and make use the recent progres s in implementing exact non-linear representations of chiral 
symmetry on the lattice 
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